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In this article, we investigate inertial modes of rigidly rotating neutron stars, i.e. modes for which 
the Coriolis force is dominant. This is done using the assumption of a fixed spacetime (Cowling 
approximation). We present frequencies and eigenfunctions for a sequence of stars with a poly- 
tropic equation of state, covering a broad range of rotation rates. The modes were obtained with 
a nonlinear general relativistic hydrodynamic evolution code. We further show that the eigenequa- 
tions for the oscillation modes can be written in a particularly simple form for the case of arbitrary 
fast, but rigid rotation. Using these equations, we investigate some general characteristics of in- 
ertial modes, which are then compared to the numerically obtained eigenfunctions. In particular, 
we derive a rough analytical estimate for the frequency as a function of the number of nodes of 
the eigenfunction, and find that a similar empirical relation matches the numerical results with 
unexpected accuracy. We investigate the slow rotation limit of the eigenequations, obtaining two 
different sets of equations describing pressure and inertial modes. For the numerical computations 
we only considered axisymmetric modes, while the analytic part also covers nonaxisymmetric modes. 
The eigenfunctions suggest that the classification of inertial modes by the quantum numbers of the 
leading term of a spherical harmonic decomposition is artificial in the sense that the largest term 
is not strongly dominant, even in the slow rotation limit. The reason for the different structure of 
pressure and inertial modes is that the Coriolis force remains important in the slow rotation limit 
only for inertial modes. Accordingly, the scalar eigenequation we obtain in that limit is spherically 
symmetric for pressure modes, but not for inertial modes. 

PACS numbers: 04.40.Dg 



I. INTRODUCTION 

Neutron star oscillations have become an interesting 
field of research because they are possible sources of grav- 
itational waves, which might be directly detected in the 
near future. Also, some are believed to undergo rota- 
tional instabilities where energy is transferred from the 
rotation of the star to the oscillation and gravitational 
radiation, which also carries away angular momentum. 
Such mechanism would not only be a strong source of 
gravitational waves, but also a possible explanation for 
the observed limitation of the rotation rate of neutron 
stars. Should gravitational waves ever be detected, we 
would gain some insight into the behavior of cold, dense 
matter, see pQ. For this, one needs accurate theoretical 
models. 

There are different families of oscillations: Spacetime 
modes (w-modes) are oscillations of the spacetime itself 
with weak coupling to the matter, see [2], and only ex- 
ist in the framework of general relativity. For pressure 
modes (f- and p-modes), the restoring forces are mainly 
pressure and gravitation. They exist in all stars, regard- 
less of entropy gradient and rotation rate. If the star 
has a radial specific entropy gradient, there also exist 
g-modes, where the restoring force is buoyancy. For iner- 
tial modes, the Coriolis force, gravitation, and pressure 
are equally important restoring forces. We will only con- 
sider the case of isentropic stellar models, for which no 
g-modes exist; according to [3J inertial modes of non- 
nonbarotropic stars behave qualitatively different. For- 
mally, inertial modes only exist in rotating stars, but for 
arbitrary low rotation rates. In the nonrotating limit, 



their frequency goes to zero and they become stationary 
currents. 

So far, astrophysical interest has been focused on a 
subclass of nonaxisymmetric inertial modes called r- 
modes, because in the absence of viscosity and nonlin- 
ear mode coupling effects, they can undergo rotational- 
gravitational instabilities at any rotation rate, see [1]. 
Recent Newtonian studies [SJ E] argue that the insta- 
bilities are in fact strongly suppressed by mode coupling. 
However, for the axisymmetric modes we extracted nu- 
merically, no rotational instabilities exist. 

Our knowledge on neutron star oscillations in general 
relativity is mainly based on perturbative studies. In [S] 
El [101 [TTJ [121 H3J [III US] oscillations of nonrotating stars 
are investigated in full relativity. Pulsations of slowly 
rotating stars are studied in [3J OH [TO HH EE! W\- In 
|211 122j . nonaxisymmetric pressure modes and r-modes 
are studied for arbitrary rotation rates, but using the 
Cowling approximation 

Most results are limited to the slow rotation approxi- 
mation, because in that limit it is feasible to decompose 
the perturbed quantities into spherical harmonics, reduc- 
ing the equations to one spatial dimension. Linearized 
evolution equations in full relativity using this formalism 
can be found in [TBI EE] [23]. In [T7], inertial modes are 
investigated using the additional assumption of a fixed 
spacetime. 

As shown in the above references, rotation couples 
spherical harmonic contributions of different quantum 
number I. This yields an infinite system of equations, 
which is truncated to I < l m to obtain numerical solu- 
tions. For the case of inertial modes however, a problem 
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arises: as detailed in [17J, the resulting equations ad- 
mit solutions only outside certain frequency bands, which 
change their location and number with the chosen trun- 
cation l m in an apparently non-converging fashion. 

In the analytic part of our work, we try to investigate 
the inertial mode problem from a different angle and di- 
rectly study the two dimensional perturbed eigenequa- 
tions. Additionally, we do not restrict ourselves to the 
slow rotation limit from the beginning, although it is used 
for some of our conclusions. We do however rely on the 
assumption of rigid rotation. 

Oscillations of neutron stars have also been studied 
using nonlinear relativistic hydrocodes, which are not re- 
stricted to the slow rotation approximation. In [24} 125] . 
pressure modes frequencies are computed within the 
Cowling approximation, and in |26j using the conformal 
flatness approximation. The latter also contains some in- 
ertial mode frequencies, but without identification of the 
eigenfunction. In [27 , a simulation of one inertial mode 
with a high amplitude in the nonlinear regime has been 
carried out in full relativity. To our knowledge, there are 
no inertial mode eigenfunctions of rapidly rotating stars 
available in the literature so far. 

This article is organized as follows: In Sec. |TT] we 
present the eigenequations and investigate the general 
characteristics of inertial modes, particularly in the slow 
rotation limit, which is derived in Sec. |IIB| In Sec. |III| 
we present numerically extracted eigenfunctions and fre- 
quencies, which are then compared to the analytic expec- 
tations. 



II. ANALYTIC PROPERTIES 

In this section, we present the eigenequations in Cowl- 
ing approximation for the case of rigid rotation. Al- 
though there is no known analytic solution, we deduce 
some general properties of inertial modes. They will be 



used in Sec. Ill to check and explain the numerical re- 
sults. 

We use the following notation: p is the rest mass den- 
sity in the fluid restframe, e the specific internal energy, 
excluding restmass, h = 1 + e + P/p the specific relativis- 
tic enthalpy, c s the speed of sound, the 4-velocity of 
the fluid, v l the 3- velocity measured by an observer with 
worldline normal to the hypersurface of constant coor- 
dinate time, W the corresponding Lorentz factor, and 
w 4 = u l /u° the advection speed of the fluid with re- 
spect to the coordinates. The 3-metric is denoted by 
gij, its volume element by ^/t~ = -^det(gij), and the 
Lapse function by a. We further use the conserved mass 
density D = y/^Wp, the conserved momentum density 
Si = y/^yW 2 phvi, and the specific angular momentum 
lip = S v /D. The angular velocity of the star measured 
by a distant observer is denoted by fl. 

We use a cylindrical type coordinate system (t, d, z, ip) 
in which the stationarity of the spacetime and the ax- 
isymmetry are manifest, i.e. all background quantities 



are invariant under translation in t- and (^-direction. We 
further require the coordinate system to be corotating, 
i.e. w l = 0. This is always possible for rigid rotation, but 
otherwise not. Spacelike indices are denoted by letters 
i,j G {d, z,(p}, while indices l,k € {d, z} denote compo- 
nents in the meridional plane. In perturbation equations, 
the perturbed quantities are generally prefixed with a 
S, while everything else refers to the unperturbed back- 
ground model. If not noted otherwise, geometric units 
where G = c = 1 are used. 



A. Linearized equations 

To obtain perturbed evolution equations, we lin- 
earize the hydrodynamic evolution equations on arbitrary 
curved spacetimes in the formulation given in [28 , which 
read 



d t D = -d s (DvP) , 

d t Si = -dj \SiVp +a^f(p- P) 81 



IQ 



(1) 
(2) 



The second equation is only valid at the arbitrarily cho- 
sen point Q, for which the function P is defined by 



P(Q) = P(Q), 



w 



h(P) — const, 



(3) 



where W is the Lorentz factor that belongs to a constant 
advection speed w ! = w t (Q). In this formulation, the ge- 
ometrical source terms which are usually written in terms 
of derivatives of the metric tensor, are combined in the 
derivative of P. For the following, it is not necessary to 
understand this formulation, which is explained in de- 
tail in [28]. The energy equation also derived there is 
redundant due to the assumption of adiabatic evolution. 

The above equations are derived assuming the stress 
energy tensor of an ideal fluid and the conservation of 
restmass. Here we further assume an isentropic stellar 
model and adiabatic evolution, both satisfying the same 
one-parametric equation of state (EOS). For the unper- 
turbed stellar model, we generally assume axisymmetry 
and rigid rotation. As shown in |28j . any stationary, isen- 
tropic, and rigidly rotating stellar model has to satisfy 



W 



const. 



(4) 



Using this as well as the axisymmetry of the unperturbed 
background model, a straightforward, but lengthy calcu- 
lation leads to the evolution equations for the perturba- 
tion 

d t SD = ~d l (D8w l ) , (5) 
8 t SS v = -ft (S v Sw l ) - kD8 v Ss, (6) 
d t SSi = D(-k8i5s + Swiftly), (7) 

where Ss = 61n(h). Note the above equations are valid 
only in corotating coordinates. The system is closed by 
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the algebraic relations 



5w h 



5s 



SSi 


= Dh 


W 

— 




a 


ss v 


= 


(( 


SD 


= °( 





w 2 



(8) 



(10) 



Eq. (|5| describes the conservation of rest mass. Because 
of the axisymmetry of the background model, there is a 
conserved angular momentum S v , described by Eq. 
To understand Eq. Q, we compute its Newtonian limit, 
obtaining 



d 5v d = 
d Sv z = 



-d d 5s 
-d z Ss, 



(ii) 

(12) 



The term 2£l6v v d, which originated from the term 
Sw^dil^ in Eq. ([7]), is nothing but the Coriolis force com- 
ponent acting in d-direction, caused by the change of ve- 
locity in 93-direction. The forces due to the perturbed 
balance of pressure and gravitational force are all con- 
tained in the term —diSs. Note the latter is only possible 
when using Eq. @. 

We now assume harmonic time dependence for all per- 
turbations SZ. Because of the axisymmetric background, 
we also assume harmonic dependence on <p, writing 



8Z(d, z,(p,t) = SZ(d, z )e i{uJt+m ^ . 



(13) 



Note to is the frequency with respect to coordinate time 
in the corotating frame; the frequency in a nonrotating 
frame is shifted by mf2. After some straightforward com- 
putations, we obtain the eigensystem 



8L 



5D 



5x k 



m—5s, 

LO 



Sx dil v - 

di (DSx 1 ) +mD-Sw v , 
x ' to 

W 



(14) 
(15) 



A k i - di5s 



2qn\v\—b l ds. (16) 
a 



Here, Sx k = Sw k /(iuj) is the fluid displacement vector, 
and 



kl — 


9ki ohh, 


(17) 


bi = 


a M » 1 ; 
01 In lu, 


(18) 




v = 


LO 

20' 


(19) 


q = 


K 

1 + m— . 

UJl v 


(20) 



The vector b l depends weakly on the rotation rate; in 
the Newtonian limit, we obtain b l — (1,0). The above 
system of equations has real coefficients. The physical 
(real- valued) solutions for the quantities Ss, 5w v , and Sx k 



therefore oscillate in phase, which is spatially constant. 
The velocity components 8w k in the meridional plane are 
phase-shifted by tt/2 with respect to 8s. 



If v 2 7^ b l bi, A k can be inverted and Eq. (16 1 can be 
solved for 5x l . Together with Eq. ( 15 ) , we obtain a second 



order scalar master equation for Ss. 



= u 2 Di]5s + di 



a 

W 2 



DC lk d k Ss 



where 



g kl - ^rb k b l , 

V 1 



b% 



(21) 

(22) 
(23) 
(24) 



A very similar formulation, valid also for differential rota- 
tion, is given in |29j . The principal part of the equation is 
determined by the tensor C k , which has the eigenvalues 
{1,— /i 2 }. One can therefore expect qualitatively differ- 
ent families of solutions for the cases v > v c and v < v c . 
From the Newtonian limit of Eq. (14 1 and Eq. (16), it 



is easy to show that for axisymmetric modes, the ra- 
tio between the d-components of the Coriolis force and 
the total net force acting on a fluid element is given by 
\jv 2 . Identifying the relativistic terms that reduce to 
the Coriolis force in the Newtonian limit, we obtain in 
the relativistic case that for v < v c , the Coriolis force 
is not only the dominant restoring force, it even has to 
be counterbalanced by the pressure force. In d-direction, 
pressure forces and fluid acceleration then become an- 
tiproportional. If v 3> v c on the other hand, the Coriolis 
force becomes unimportant. This is exactly the criterion 
that distinguishes inertial modes and pressure modes. 



B. Slow rotation limit 

In the following, we investigate the slow rotation limit 
O — s- of the eigenequations. It is important to note that 
for inertial modes, the oscillation frequency w is propor- 
tional to the rotation rate O. Therefore, just ignoring 
terms of higher order in £1 in the eigenequations would 
be wrong. First, we have to replace to by 2£lv. Then we 
take the limit f2 — > Q,v — > i/q, where isq is some finite 
value. 

For pressure modes on the other hand, the oscillation 
frequency to approaches a finite value loq for Q — > Q. 
Hence the proper slow rotation limit is given by f2 — > 
0, lo —>■ loq > 0, v — * 00. 

In this case, C kl — > g kl . 
limit then becomes 



Eq. (21 1 in the slow rotation 



,2 " \ 

g<pip / 



(25) 
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Using Eq. (13) and the axisymmetry of the background, 



this can also be written in the covariant form 



= Ldn-^Ss 



(Vj (a 2 p)) \7 l Ss + a 2 pASs. 



(26) 



Neglecting rotational corrections of second order in the 
rotation rate, the quantities p, c s , and a are given by the 
spherically symmetric TOV solution, see [301 l3"T] . The 
above equation becomes invariant under rotation and can 
be solved by the ansatz 



5s(r,9,<p)=5s(r)Yr(e,v). 



(27) 



It has been shown that there is an infinite discrete set 
of solutions for given I, \m\ < I, exactly one for a given 
number n of radial nodes. The frequency grows without 
bound with increasing n or I. 

For the incrtial modes, such an ansatz is not possible. 
The reason is that in the slow rotation limit for inertial 
modes, we obtain 



C 



kl 



/~ikl Kl 

°o — 9 ~ 



kl 



1 



;b k b l 



(28) 



The second term is finite, which means that other than 
the star profile, the equations themselves are not even 
approximately spherically symmetric. This is reflecting 
the fact that the Coriolis force, which is a directional 
force, remains dominant for arbitrary slow rotation rates. 

In perturbative mode calculations like [T7], the scalar 
quantities are usually decomposed into spherical har- 
monic contributions fn( r )Yi m (®)- For pressure modes, 
this is a natural choice since in the nonrotating limit, one 
such term is enough to describe the solution, while for 
moderate rotation rates, the others remain small correc- 
tions. 

For inertial modes on the other hand, one cannot a 
priori expect that the above decomposition is more nat- 
ural than any other expansion, e.g. as twodimensional 
Fourier series. As will be shown in Sec. |III| the angular 
dependency of axisymmetric inertial mode eigenfunctions 
does indeed not bear resemblance to spherical harmonics. 
This is not to say that spherical harmonic decomposition 
is a bad choice, but to stress that inertial modes have a 
completely different structure than pressure modes. 

The cigencquations we obtain for inertial modes in the 
slow rotation limit are given by 





V 



1 



fjm 2 5s+-di(Da 2 C$ l d k 5s), 



2 ? 

Mo" 



(29) 
(30) 



In particular for m = 0, the equation is quite simple and 
might be suitable for direct numerical solution. Note 
however that the unknown v enters the equations in a 
nonlinear manner as a downside. 



From Eq. ( 14 1 to Eq. ( 16 1, it follows that in the incrtial 



mode slow rotation limit, 

o(Ss) = o(nsw l ) = o(n5w v ) = o(n 2 Sx i ). (31) 



The density perturbation thus vanishes for a finite ve- 
locity amplitude in the slow rotation limit. The physi- 
cal interpretation is that pressure and Coriolis force are 
proportional for inertial modes, and the latter is propor- 
tional to the rotation rate. The mass current has to be- 
come divergence-free. For axisymmetric oscillations, our 
numerical results show that this is realized by fluid mo- 
tions which correspond to a sum of convection-like mo- 
tions in the meridional plane and differential rotation, 
both of course with a harmonic time dependence and 
phase shifted by ir/2. 

Note that for an arbitrary small but fixed velocity am- 
plitude, the fluid displacement vector 8x % goes to infinity 
in the slow rotation limit. However this does not mean 
the linear approximation is not valid anymore to describe 
such an oscillation, since the Eulerian perturbations re- 
main all finite. Only for the description of the movement 
of a fluid element one would have to integrate the veloc- 
ity field along its path, which will result in the aforemen- 
tioned motions. For slow rotation rates, one can thus 
have arbitrary many turnarounds during one oscillation 
period. 

We now briefly discuss what happens in the presence 
of a small, but finite entropy gradient. Lets assume the 
fluid moves adiabatically with exactly the same velocity 
field as in the isentropic case. The density perturbation 
would then also be the same. The Eulerian perturbation 
of specific entropy would be proportional to the fluid dis- 
placement vector, if the latter is small, but remain finite 
in any case, since the fluid elements only move around 
inside the star as argued before. Due to the entropy per- 
turbation, there is an additional pressure perturbation 
compared to the isentropic case. The size of this pertur- 
bation has to be compared to the pressure perturbation 
corresponding to the one of the isentropic inertial mode 
oscillation. If it is much smaller, it can be regarded as 
a small correction to the isentropic inertial mode oscil- 
lation, which is still an approximate solution. If it is 
much bigger on the other hand, the oscillation must have 
a completely different structure than in the isentropic 
case. 

If the entropy gradient is small enough, we thus expect 
three different regimes for finite amplitude incrtial mode 
oscillations: a small amplitude regime where the struc- 
ture of the modes changes drastically, a medium ampli- 
tude regime where the modes are similar to the isentropic 
case, and the nonlinear regime. The first case is captured 
by linear perturbation theory, where all oscillations are 
arbitrary small by definition. Indeed, [3] found qualita- 
tively different solutions in the case of an entropy gradi- 
ent. Wether there is a medium regime before the non- 
linear regime starts depends on the size of the entropy 
gradient as well as on the rotation rate; demanding that 
the velocity perturbation is bound by some finite value, 
the maximum pressure perturbation of isentropic incr- 
tial modes is proportional to the rotation rate. Thus we 
expect that finite amplitude inertial mode oscillations of 
rapidly rotating stars are more robust to entropy gradi- 
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ents than in slowly rotating stars. 



C. Boundary conditions 

We now investigate the behavior of the cigenequations 
at the surface. For this, we assume that Ss and its first 
and second derivatives remain finite at the surface. This 
assumption is compatible with our numerical results. It 
is easy to show that if there is a surface at all, \d l p\ S> p 
near the surface. Neglecting the terms directly propor- 
tional to p in Eq. (21), we arrive at 



ujHsV = C kl (d lP )(d k 5s), 



V 



w 

a 



2n^ 



qb l dip 



(32) 
(33) 



For the following, we assume that p/c 2 remains finite near 
the surface, as it is the case for the polytropic EOS used 
in our numerical computations. For the axisymmetric 
case where q = 1, the slow rotation limit for inertial 
modes then leads to the condition 



if)«d k 8s = 0, V* = C diP 



(34) 



at the surface. 

To investigate the behavior of if) , we decompose dip 
into the normalized eigenvectors of C l k , which are given 



by 



b l /\b\, with eigenvalue — p 2 , and e l ± defined by 



e\bi = 0, | ex | = 1, corresponding to the eigenvalue +1. 
The eigenvectors ey and e± are roughly (in the Newto- 
nian limit exactly) orthogonal, respectively parallel, to 
the rotation axis. For small v and hence small p, ip l 
will point in the direction of e-j away from the axis, ex- 
cept near the equatorial plane, where e\dip rs 0. In that 
case, nodes crossing the surface at some distance from 
the equatorial plane will be approximately parallel to the 
rotation axis. 



D. Small wavelength limit 



Since the numerically extracted inertial modes shown 



in Sec. Ill expose increasingly regular patterns with grow- 
ing number of nodes, it is instructive to investigate the 
limit of small wavelengths. Let's assume there is a so- 
lution for which the first and second derivatives of Ss, 
normalized to its maximum amplitude, are much greater 
than the normalized derivatives of the background quan- 
tities. Treating the background quantities as constant in 
the neighborhood of some point inside the star, Eq. (21 ) 
becomes 







USs + C kl d k diSs, 



U = uo z 



,w 



This can be solved by the ansatz 

Ss = 8s e lk ' Xl . 



(35) 
(36) 

(37) 



We now de comp ose k l into normalized eigenvectors of 
C l k , see Sec. II C writing k l — fcjej + k±e l ± . It follows a 
dispersion relation 



p- 



k 2 ± -U 



k 2 r 



(38) 



Globally, the eigenfunctions cannot be described as plane 
waves. Nevertheless, we can roughly estimate the lo- 
cal second derivatives in each direction from the number 
of nodes along the (entire) rotation axis and equatorial 
plane, setting fcj_ « n r ir/(2Rp) and k-j ~ n e ir/R e , where 
R p and R e are the polar and equatorial coordinate radius 
of the star. 

For the axisymmetric case, we can also neglect U in 
the slow rotation approximation for inertial modes, which 
yields the relation 



R e n r 
2R V n e 



(39) 



We stress that the above approximation is quantitatively 
very bad. For example, v is a global constant while v c 
varies by 30 % in the stellar models we investigated. 

However, qualitatively we expect that frequencies of 
inertial modes increase with the number of nodes along 
the rotation axis, but decrease with the number of nodes 
along the equatorial plane. We further expect that for 
high order modes v is bounded by some i7 ci which is of 
similar magnitude as the values of v c inside the star. Fi- 
nally, for small p, there should be nodes roughly parallel 
to the rotation axis. 



E. Beyond Cowling 

In the following we will briefly discuss the emission of 
gravitational waves and the quality of the Cowling ap- 
proximation, for inertial modes in slowly rotating stars. 



A consequence of Eq. ( 31 1 for estimating the gravita- 



tional waves emitted by inertial modes is that the dom- 
inant contribution in the quadrupole formula is given 
by the current terms, not the density terms. This can 
be seen as follows: the second time derivative of the 
mass quadrupole is given by I = (p) m ~ ^ 2 (Sp) m , 
where the brackets denote the integral operators given 
in |32j . while for the current quadrupole we have S = 
(df{pv 1 )) = fl 2 (v l Sp + pSv 1 ) . Since 8p ~ Ss, v ~ fl 

and O(Ss) = 0(n8w l ), we get 0(1) - 0(n 2 6s), but 
O(S) - 0(QSs). 

However, in the slow rotation limit the gravitational 
waves produced by inertial modes become negligible since 
their frequency goes to zero; for the luminosity, we obtain 
L - | S | 2 - n 6 (Sv*) 2 e , or O(L) - 0(n 6 E k ), where E k is 
the kinetic energy in the oscillation. Therefore gravita- 
tional radiation should not play a role for the dynamics 
of inertial mode oscillations of slowly rotating stars. 
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Assuming that dissipative and mode-coupling effects 
limit the velocity perturbation amplitude to \Sv l \ < vm, 

H, 6 for the gravita- 



(a) xle-4 



' M 



it follows an upper bound L m 

tional wave luminosity produced by any inertial mode. 

Besides gravitational radiation, the Cowling approxi- 
mation neglects already on the Newtonian level the per- 
turbation of the gravitational potential. It is unclear how 
this affects inertial modes. Although the density pertur- 
bation and hence the perturbation of the potential goes to 
zero in the slow rotation limit for a fixed velocity ampli- 
tude, the pressure perturbation and the Coriolis force do 
the same and cannot be neglected. For pressure modes, 
it is known that the Cowling approximation can be quite 
inaccurate, as shown in [26]. In [33], it has been argued 
that for the case of pressure modes of nonrotating stars 
the Cowling approximation becomes exact in the low fre- 
quency limit. Note however that the Cowling approxi- 
mation discussed there is a refinement where the metric 
component g rt is also perturbed; see [34 for a compari- 
son of the two variants. Anyway, these results cannot be 
directly carried over to inertial modes of rotating stars. 
We thus believe that the Cowling approximation might 
be more accurate for inertial modes than for pressure 
modes, but this can only be quantified by a fully rela- 
tivistic study. 



III. NUMERIC RESULTS 

The aim of this section is to give a qualitative overview 
how inertial modes actually look like, to investigate the 
influence of the rotation rate, and to provide accurate 
frequencies for reference in future works. 



A. Method 

To obtain frequencies and eigenfunctions, we perturb 
an equilibrium stellar model with a small trial pertur- 
bation and then evolve the general nonlinear hydrody- 
namic equations in time, using the PIZZA code described 
in [28 : 35 . We use the Cowling approximation, i.e. the 
spacetime is kept fixed. This saves a lot of resources, 
since we only need to evolve the interior of the star, but 
not a huge volume of surrounding spacetime. The am- 
plitudes of the perturbations are generally in the linear 
regime, with velocities in the range 10~ 7 to 10~ 4 c. We 
use the nonlinear pizza code instead of a linear code only 
because it is available and well tested. 

From the Fourier spectra of the density and velocity 
variations at some sample point inside the star, we deter- 
mine the frequencies of the excited oscillations. We then 
perform a Fourier analysis for one of those frequencies at 
every point inside the star, to obtain a first approxima- 
tion to the eigenfunctions. The real part of the complex 
result is then used as initial perturbation in a second sim- 
ulation. This step is repeated until other modes are suffi- 
ciently suppressed, like in the time evolution and spectra 
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(b) 




FIG. 1: Time evolution (a) and Fourier spectrum (b) of ve- 
locity in ^-direction at r — 3 km, 6 — n/2, for a simulation of 
stellar model BU3, which was perturbed using the eigenfunc- 
tion of inertial mode 122 • The resolution of the simulation is 
150 x 150 points. The damping is purely numerical. 



shown in Fig. [T] This method called mode recycling has 
already been used in |26j . 

As shown in Sec. |II A| the velocity components in the 
meridional plane are phase shifted by 7r/2 against the per- 
turbations of specific energy and velocity in ^-direction. 
Apart from the global complex phase, the eigenfunctions 
are purely real. To obtain the eigenfunctions from the 
complex-valued result of the Fourier analysis, we first 
compute a weighted mean phase by 



I %/7P( ar g(^ £ c) mod 7r) d 2 x 



(40) 



where 5e c is the result of the Fourier analysis for the 
specific energy perturbation. The eigenfunctions are then 
obtained from Se = R(e-^°5e c ), Sw* = ^(e'^Sw^), 
and Sw l = Sy(e-^°5u£). 

We also compute the quantity cos(arg((5e c ) — 0o), which 
should result in ±1, with a sign flip at the nodes of 
the eigenfunction Se. This is used to estimate the accu- 
racy of the extracted mode, in particular to determine if 
more mode recycling steps are necessary. Examples how 
this measure looks like in our computations are shown in 
Fig. [I 
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FIG. 2: Uniformity of the global complex phase under ideal 
circumstances (a) and for a problematic case (b). Color coded 
is the quantity cos(arg(5e c ) — cj>o). 



The method described above has been successfully ap- 
plied to pressure modes in [28, 35 . For inertial modes, 
the scheme is also usable, but computationally less ef- 
ficient. The reason is that the ability of the Fourier 
analysis to separate different frequencies and the cor- 
responding eigenfunctions depends on the time interval 
given by the evolution time. For pressure modes, this 
is not a problem because their frequencies are well sep- 
arated. For inertial modes, there seem to be infinitely 
many modes inside a given frequency interval. Although 
most of them are of high order and cannot be resolved 
by simulations with a finite resolution, the frequencies 
of resolvable modes can still become much closer than 
those of the pressure modes, which makes it more diffi- 
cult to distinguish them by Fourier analysis. Note that 
high order inertial modes can have similar frequencies as 
low order inertial modes, in contrast to pressure modes. 

The direct approach is to use very long evolution times. 
Since the usable evolution time is effectively limited by 
the timescale of the numerical damping, one would also 



need higher resolutions. The brute force method is not 
only computationally inefficient, it is also not guaran- 
teed to work, since with increasing resolution more modes 
with similar frequencies are becoming resolvable. 

However, there is a simpler method to compute at least 
the lower order modes, by using the numerical errors 
causing the damping of the oscillations to our advan- 
tage. Higher order modes experience a stronger numer- 
ical damping, and are therefore suppressed by a certain 
factor during every mode recycling iteration, no matter 
how small the difference in frequency to the low order 
mode we want to extract might be. This method to sep- 
arate modes works best at medium resolutions around 
150-200 points per stellar radius; for higher resolutions, 
the numerical damping becomes too weak, while for low 
resolutions, the accuracy of the results is not sufficient. 

In practice, we use up to ten mode recycling steps, evo- 
lution times between 20-300 ms, and resolutions between 
150-300 points per stellar radius for inertial modes. For 
pressure modes, one to three iterations, evolution times 
around 10-20 ms, and resolutions around 100 points were 
usually sufficient. 

The aforementioned complications also make it hard 
to properly compute the error of the inertial mode eigen- 
functions. The error now depends not only on the accu- 
racy of the simulation, which is investigated in [2"51 135] . 
but also on the extent to which modes with frequencies 
in the same Fourier bin could be suppressed. This can- 
not be estimated from the Fourier spectra. However, if 
the eigenfunction of a neighboring mode present in the 
simulation significantly differs from the mode of interest, 
it would show up as an error in the uniformity of the 
global complex phase. We cannot separate modes with 
very similar frequencies and eigenfunctions, however. 

Based on comparisons of different simulations of the 
same mode, we estimate the error of the eigenfunctions 
to 5 %. Our estimate for the error of the frequency is 
Sf = 0.01/ + AF, where AF is the width of the Fourier 
bin corresponding to the evolution time. 



B. Inertial mode eigenfunctions 

Using the method described in the previous section, 
we computed the frequencies and eigenfunctions of five 
low order axisymmetric inertial modes for a sequence of 
rigidly rotating stars with various rotation rates and a 
fixed central density of 7.9056 x 10 17 kg/m 3 . The equa- 
tion of state is a polytrope 



Pp \PpJ 

with r = 2, p p = 6.1760 x 10 18 kg/m 3 . Our stellar 
models are named BU2-BU4 and BU6 in [2! [26], where 
the corresponding pressure modes are investigated. Ad- 
ditionally, we use a model "BUS" with a slow rotation 
rate of 246 Hz. The fastest rotating model has a rotation 
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rate of 792 Hz, and an axes ratio of 0.7. For comparison: 
the mass-shedding limit for the given EOS and central 
density is reached for an axis ratio of 0.58. Deviations 
from the slow rotation limit should therefore become vis- 
ible. We restricted ourselves to equatorially symmetric 
simulations, in order to save computational resources. 

Our results for the axisymmetric eigenfunctions of 
model BU3, which has a rotation rate of 590 Hz and an 
axis ratio of 0.85, are shown in Fig. [3] to Fig. [7] As pre- 
dicted in Sec. |II A| inertial modes have a completely dif- 
ferent structure than pressure modes, at least for the ax- 
isymmetric case. For comparison, pressure mode eigen- 
functions for model BU3 can be found in [35]. The de- 
pendency of the patterns on the rotation rate is weak for 
all modes. Fig. [8] shows one of them at two different ro- 
tation rates. Obviously, the angular dependency is not 
given by a spherical harmonic plus rotational corrections 
proportional to f2 2 , like it is the case for pressure modes. 

For practical purposes, we classify the axisymmetric 
inertial modes by the number of nodes of the eigenfunc- 
tion 5s along the entire rotation axis and the equatorial 
plane, e.g. 123 has two nodes along the (entire) axis, and 
three along the equatorial plane. It is unclear wether the 
number of nodes in each direction uniquely determines 
the mode; for the modes we extracted so far, this is the 
case. 

For some simulations, the global complex phase re- 
mained a little fuzzy in the regions where two nodes come 
close, as shown in Fig. [2] This could easily be explained 
by the fact that the phase is much more sensitive to errors 
in regions where the eigcnfunction is nearly zero. How- 
ever we cannot rule out the possibility that our modes 
are really superpositions of several modes with very sim- 
ilar frequencies whose eigenfunctions differ significantly 
only in the aforementioned regions. 

Although we performed the mode recycling steps re- 
quired to obtain a clean eigenfunction only for modes 
with two nodes along the axis, due to limited computa- 
tional resources, we also encountered patterns with more 
nodes along the axis. 

The nodes of the extrated inertial mode eigenfunctions 
tend to be oriented roughly parallel and orthogonal to 
the rotation axis, which becomes more pronounced with 
increasing number of nodes crossing the equatorial plane. 
This corresponds well to the analytical prediction from 
Sec. |HDj 



(a) 
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In order to check prediction Eq. ( 34 1 for the direction 
of the gradient of Ss in the slow rotation limit, we com- 
puted ijr using central finite differences from the same 
background model data used as initial data in the simu- 
lations. The result is compared to the numerical eigen- 
functions in Fig.|3]to Fig. [7] As one can see, the analytic 
prediction nicely matches the numerical results. 

Another visible difference of the extracted inertial 
modes to pressure modes is that the velocity fields are 
vortexlike and tangential to the surface. This reflects the 
fact that in the slow rotation approximation, the den- 
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FIG. 3: Eigenfunction of inertial mode 121, computed for stel- 
lar model BU3. The cuts show the upper half of the merid- 
ional plane, (a) Eigenfunction Ss. With the exception of the 
star surface, the lines represent equidistant isocontours. Neg- 
ative values correspond to dashed lines, and the thick solid 
line marks the node of the eigenfunction. The shaded area 
represents amplitudes below 5% of the maximum one, which 
is our guess for the accuracy of the node location. The arrows 
at the surface represent the analytic prediction for the direc- 
tion of the isopotential lines of Ss in the slow rotation approxi- 
mation, (b) Eigenfunction of momentum density components 
S l in the meridional plane. The solid lines are the nodes of 
Ss. 



sity change vanishes, as shown in Sec. |IIB| and the mass 
current must become divergence-free. 

As a consequence, the oscillation amplitude Sx l of the 
fluid displacement vector can become arbitrary large in 
the slow rotation limit without leaving the linear regime. 
The fluid then performs convection-like motions which 
reverse periodically, without changing the star profile. 
Note that for a fixed velocity perturbation amplitude, 
the fluid displacement vector 5x l has to go to infinity in 
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FIG. 4: Like Fig. [3] but for the i 22 mode. 



FIG. 5: Like Fig. [3] but for the i 2 3 mode. 



the slow rotation limit due to Eq. (31 1 



C. Inertial mode frequencies 

Together with the cigenfunctions, we extracted the fre- 
quencies. In most cases, other modes were suppressed 
strongly enough to allow a direct fit of an exponentially 
damped oscillation to the time evolution of the velocity 
components at some sample point. Otherwise, the fre- 
quencies were extracted from the corresponding Fourier 
spectrum. 

The frequencies are given in Tab. [I] As shown in 
Sec. |II A| the frequencies of inertial modes should be pro- 
portional to the rotation rate at leading order. The ac- 
tual ratio between the two depicted in Fig. [9] is indeed 
compatible with a correction term to v of second order 
in the rotation rate. By fitting the data accordingly, we 
obtain extrapolated values for the slow rotation limit, 
which is inaccessible to our numerical method. The re- 



sulting values are given in Tab.|TTJ Note that our sequence 
of stellar models does not correspond to the same star at 
different rotation rates since we fixed the central density 
instead of the total baryonic mass. However, the leading 
correction term to the total mass is also of second order. 

The modes considered here satisfy an interesting em- 
pirical relation: for a given stellar model, there is a nearly 
linear relation between the quantity 1/fi introduced in 
Sec. II A and the number of nodes along the equatorial 
plane. Note that n is not constant inside the star, and the 
accuracy of the linear parametrization slightly depends 
on the choice of the position where /x is computed (which 
has of course to be the same for all modes) . For model 
BU3, Fig.[lO]shows the minimum and maximum values of 
/i -1 as well as the Newtonian value we define by setting 
v c = 1. The latter can be conveniently computed from 
the rotation rate and frequency alone, without knowl- 
edge of the stellar model. The relation also holds for 
the other stellar models and for the extrapolated values 
in the slowly rotating case, but with a slightly different 
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FIG. 6: Like Fig. [3] but for the 124 mode. 
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FIG. 7: Like Fig. [3] but for the i 25 mode. 



slope and offsets. 

The analytical estimate Eq. ( 39 1 for higher order 
modes, see Sec. |IID| which is also shown in Fig. [10] 
matches the data only by order of magnitude. This is not 
surprising considering the drastic approximations used in 
the derivation. However, the fact that the data fulfills an- 
other linear relation quite accurately could be a hint for 
the existence of a better analytical estimate. 

As a further check, the 121 inertial mode has been ex- 
tracted by [3E], using the COCONUT code described in 
[261 137] . Since only a single mode recycling iteration was 
used, we estimate the errors of the eigenfunctions around 
20 %, judged from the observed phase errors and our 
previous mode recycling experiences with inertial modes. 
However, the eigenfunctions expose the same qualitative 
structure shown in Fig. [3] allowing us to identify the 
mode. The frequency of 634 Hz matches our result within 
2 %. 
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FIG. 8: Comparison of the eigenfunction of inertial mode 122 
at rotation rates 246 Hz (a) and 792 Hz (b). See Fig. [3] for a 
description of the plot. 
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FIG. 9: Ratio v = lo/2Q versus rotation rate Fr for a 
sequence of stellar models with fixed central density. The 
shaded region marks the range of v c inside the star. The 
solid lines are fits of the function vq + U2Q 2 to the data. The 
coefficients are given in Tab. [Hi 



Star model 


Fr/Uz 


Mode 


//Hz 


Sf/f 


V 


BUS 


246 


z'21 


252.3 


2.0 % 


0.5127 






122 


195.2 


2.8 % 


0.3967 


BU2 


487 


%2\ 


507.7 


6.0 % 


0.5209 






222 


399.8 


1.9 % 


0.4102 






223 


324.4 


4.9 % 


0.3329 






2*24 


272.1 


4.7 % 


0.2792 






«25 


234.2 


4.6% 


0.2403 


BU3 


590 


221 


621.8 


1.6 % 


0.5262 






222 


496.4 


2.7 % 


0.4201 






223 


408.6 


7.2 % 


0.3458 






224 


342.9 


4.7 % 


0.2902 






*25 


296.2 


6.7% 


0.2507 


BU4 


673 


221 


717.0 


2.6 % 


0.5325 






222 


581.0 


3.9 % 


0.4315 






223 


481.5 


6.2 % 


0.3576 






224 


410.4 


7.1 % 


0.3048 






«25 


356.6 


5.7% 


0.2649 


BU6 


792 


221 


865.3 


6.8 % 


0.5462 






*22 


728.2 


5.6 % 


0.4597 






*23 


621.2 


3.1 % 


0.3921 






224 


535.5 


5.7% 


0.3380 






125 


469.0 


4.6 % 


0.2960 



TABLE I: Frequencies / of various inertial modes for the 
sequence of stellar models described in the main text. 5f is 
our estimate for the error of /. Fr denotes the rotation rate 
of the star. 



Mode 


221 


«22 


«23 


«24 


«25 




0.508 


0.386 


0.294 


0.239 


0.203 




0.058 


0.109 


0.152 


0.153 


0.145 



TABLE II: Coefficients of the parametrization v = + 
obtained by a fit to the numerical results for each inertial 
mode. 



O O ff/ 

— Linear fit 
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2 3 4 

Number of nodes 

FIG. 10: Dependency of 1/ fi on the number of nodes crossing 
the equatorial plane, for the inertial modes of stellar model 
BU3 with exactly two nodes crossing the (entire) rotation 
axis. Plotted is the value in the Newtonian limit, where \x 
is computed setting v c = 1, as well as the minimum and 
maximum values fio, A*i in the actual stellar model. The solid 
line is a linear fit to the Newtonian values, the dotted line is 
the estimate Eq. (391 computed with v c = 1. 
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IV. SUMMARY 

In this work, we extracted the frequencies and eigen- 
functions of several axisymmetric inertial modes of 
rigidly rotating neutron star models with a polytropic 
equation of state, for slow to medium rotation rates. 
For this, we used a nonlinear code for the evolution of 
ideal fluids in arbitrary spacetimes in conjunction with 
the mode recycling technique for the extraction of single 
oscillation modes. The spacetime is kept fixed, i.e. we 
use the Cowling approximation. 

The extracted frequencies are proportional to the rota- 
tion rate at leading order, and well below two times the 
rotation rate. The spectrum of mode frequencies seems 
quite dense, which made it more difficult to extract single 
inertial modes, compared to pressure modes. Although 
the slow rotation limit is not accessible to our numerical 
method, we computed the frequencies in the slow rota- 
tion limit using extrapolation. 

The scalar eigenfunctions of the axisymmetric inertial 
modes exhibit a checkerboard like structure. In contrast 
to pressure modes, it seems unnatural to classify axisym- 
metric inertial modes by the dominant term of a decom- 
position into spherical harmonics; which term is domi- 
nant depends on the radius, and also globally there is 
no strongly dominant term. That does not mean such 
decomposition is not a useful technique, but one should 
keep it in mind when talking e.g. of a I — 2, m = inertial 
mode. Note the former might not apply for nonaxisym- 
metric modes, which we did not extract. For practical 
purposes, we used the number of nodes along equato- 
rial plane and rotation axis to classify the axisymmetric 
modes. 

In the analytic part of our work, we derived a simple 
scalar eigenequation describing the eigenfunctions in the 
case of rigid rotation. Using this equation, we investi- 
gate the behavior of inertial modes at the surface, which 
agrees well with our numerical results. Further, we de- 
rived an approximate relation between the frequency of 
the mode and the number of nodes along the equatorial 
plane and the rotation axis. This relation predicts that 
the frequency grows with the number of nodes along the 
rotation axis, but decreases with the number of nodes 
along the equatorial plane. The latter indeed holds for 
the numerical results. The first part has not been vali- 
dated yet since so far we systematically investigated only 
modes with two nodes along the axis. 

The relation was derived for the case of higher order 
modes, and even then the expected accuracy is around 
30 % at most. Surprisingly, the numerical results fulfill a 
similar relation quite accurately, even for the lowest order 
modes. This leads us to speculate wether this empirical 
relation could possibly be derived analytically. 

We also investigated the slow rotation limit. Since 
inertial mode frequencies are proportional to the rota- 
tion rate, in contrast to pressure modes, we obtain two 



different eigenequations describing pressure and inertial 
modes. The difference is due to terms related to the Cori- 
olis force, which remains dominant in the slow rotation 
limit for inertial modes, but becomes negligible for pres- 
sure modes. Accordingly, the second order scalar partial 
differential eigenequation we obtain in that limit is in- 
variant under any rotation for pressure modes, but for 
inertial modes there is a preferred direction given by the 
rotation axis. 

We have shown that in the slow rotation limit, gravi- 
tational radiation of inertial modes in isentropic stars is 
mainly due to the mass currents instead of the density 
perturbations, at least on the level of the quadrupole for- 
mula. Further, the gravitational radiation should become 
negligible for the dynamics in the slow rotation limit, as 
it's luminosity goes like f2 6 for a fixed kinetic energy of 
the oscillation. 

We have no quantitative estimate how strongly inertial 
modes are affected by neglecting the perturbation of the 
gravitational potential. The available results e.g. of [25] 
for pressure modes cannot be applied to inertial modes 
since the relations between the magnitudes of different 
quantities are completely different. To estimate the ac- 
curacy of the Cowling approximation for inertial modes, 
we need fully relativistic studies. 

Since the axisymmetric inertial modes investigated 
here are not excited by rotational instabilities and are 
weak emitters for all but the most rapidly rotating stars, 
for which their frequency becomes comparable to those of 
pressure modes, a direct detection of the resulting grav- 
itational waves seems unlikely. We are also not aware of 
any other mechanism which would allow to detect those 
modes. However, given our limited knowledge of neutron 
star physics, this might quickly change. We believe it is 
important to know the complete mode spectrum of neu- 
tron stars for future studies, e.g. of mode coupling effects 
or the interaction of oscillations with the crust and the 
magnetosphere . 

Finally, we gave heuristic arguments why inertial 
modes should have a different structure in the presence of 
an entropy gradient, but only if the pressure perturbation 
of the isentropic oscillation becomes small in comparison 
to the one caused by the entropy gradient; as a conse- 
quence, oscillations with a given kinetic energy are more 
sensible to entropy gradients for slowly rotating stars. 



Acknowledgments 

This work was supported by the collaborative research 
center 'SFB-Transregio 7 Gravitationswellenastronomic" 
of the German Science Foundation (DFG) . Computations 
have been carried out mainly on the "Pioneer" cluster of 
the University of Tubingen. Thanks go to Harald Dim- 
mclmcier and Kostas Kokkotas for fruitful discussions. 



13 



[1] N. Andersson and K. D. Kokkotas, Mon. Not. R. Astron. 

Soc. 299, 1059 (1998). 
[2] K. D. Kokkotas and B. F. Schutz, Mon. Not. R. Astron. 

Soc. 255, 119 (1992). 
[3] K. H. Lockitch, N. Andersson, and J. L. Friedman, Phys. 

Rev. D. 63, 024019 (2000). 
[4] N. Andersson, Astrophys. J. 502, 708 (1998). 
[5] L.-M. Lin and W.-M. Suen, Mon. Not. R. Astron. Soc. 

370, 1295 (2006). 
[6] J. Brink, S. A. Teukolsky, and I. Wasserman, Phys. Rev. 

D. 70, 121501(R) (2004). 
[7] P. Arras, E. E. Flanagan, S. M. Morsink, A. K. Schenk, 

S. A. Teukolsky, and I. Wasserman, Astrophys. J. 591, 

1129 (2003). 

[8] K. S. Thorne and A. Campolattaro, Astrophys. J. 149, 
591 (1967). 

[9] R. Price and K. S. Thorne, Astrophys. J. 155, 163 (1969). 
[10] K. S. Thorne, Astrophys. J. 158, 1 (1969). 
[11] K. S. Thorne, Astrophys. J. 158, 997 (1969). 
[12] A. Campolattaro and K. S. Thorne, Astrophys. J. 159, 
847 (1970). 

[13] J. R. Ipser and K. S. Thorne, Astrophys. J. 181, 181 
(1973). 

[14] S. Detweiler and L. Lindblom, Astrophys. J. Suppl. Ser. 
53, 73 (1983). 

[15] S. Detweiler and L. Lindblom, Astrophys. J. 292, 12 
(1985). 

[16] J. Ruoff, A. Stavridis, and K. D. Kokkotas, Mon. Not. 

R. Astron. Soc. 332, 676 (2002). 
[17] J. Ruoff, A. Stavridis, and K. D. Kokkotas, Mon. Not. 

R. Astron. Soc. 339, 1170 (2003). 
[18] Y. Kojima, Phys. Rev. D. 46, 4289 (1992). 
[19] Y. Kojima, Astrophys. J. 414, 247 (1993). 
[20] K. H. Lockitch, J. L. Friedman, and N. Andersson, Phys. 

Rev. D. 68, 124010 (2003). 
[21] S. Yoshida and Y. Eriguchi, Astrophys. J. 490, 779 

(1997). 

[22] S. Yoshida, S. Yoshida, and Y. Eriguchi, Mon. Not. R. 

Astron. Soc. 356, 217 (2005). 
[23] A. Stavridis and K. D. Kokkotas, Int. J. Mod. Phys. B 

14, 543 (2005). 

[24] N. Stergioulas, T. A. Apostolatos, and J. A. Font, Mon. 
Not. R. Astron. Soc. 352, 1089 (2004). 

[25] J. A. Font, H. Dimmelmeier, A. Gupta, and N. Ster- 
gioulas, Mon. Not. R. Astron. Soc. 325, 1463 (2001). 

[26] H. Dimmelmeier, N. Stergioulas, and J. A. Font, Mon. 
Not. R. Astron. Soc. 368, 1609 (2006). 

[27] N. Stergioulas and J. A. Font, Phys. Rev. Lett. 86, 1148 
(2001). 

[28] W. Kastaun, Phys. Rev. D. 74, 124024 (2006). 
[29] J. R. Ipser and L. Lindblom, Astrophys. J. 389, 392 
(1992). 

[30] J. Oppenheimer and G. Volkoff, Phys. Rev. 55 (1939). 
[31] J. B. Hartle, Astrophys. J. 150, 1005 (1967). 
[32] K. S. Thorne, Rev. Mod. Phys. 52, 299 (1980). 
[33] L. S. Finn, Mon. Not. R. Astron. Soc. 232, 259 (1988). 
[34] L. Lindblom and R. J. Splinter, Astrophys. J. 348, 198 
(1990). 

[35] W. Kastaun, Ph.D. thesis, Eberhard-Karls Univer- 
sity Tubingen (2007), URL |http://tobias-lib.ub. 
uni-tuebingen . de/volltexte/2007/28037] 



[36] H. Dimmelmeier, private communication (2008). 
[37] H. Dimmelmeier, J. Novak, J. A. Font, J. M. Ibafiez, and 
E. Muller, Phys. Rev. D. 71, 064023 (2005). 



